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ABSTRACT 


This  paper  discusses  the  estimation  of  the  Doppler  frequency  of  an  incoming  RF 
signal  that  is  corrupted  by  additive  white  Gaussian  noise.  The  estimation  is  performed 
using  the  sampled  in-phase  and  quadrature  components  of  the  demodulated  signal.  The 
Cramer-Rao  lower  bounds  for  the  minimum  variance  estimators  as  well  as  the  maximum- 
likelihood  (ML)  estimator  are  derived.  The  mean-square  error  of  the  ML  estimator  is 
then  studied. 


RESUME 

Get  article  discute  de  I’estimation  de  la  frequence  Doppler  d’un  signal  RF  cor- 
rompu  par  du  bruit  blanc  Gaussien  et  additif.  L’estimation  s’est  faite  en  utilisant  les 
composantes  en  phase  et  en  quadrature  du  signal  demodule.  Les  homes  inferieures  de 
Cramer-Rao  (C-R)  pour  les  estimateurs  a  variance  minimale  ainsi  que  I’estimateur  a  prob- 
abilite  maximale  sont  derives.  L’erreur  quadratique  moyenne  de  I’estimateur  a  probabilite 
maximale  est  ensuite  etudiee. 


A0««88lon  for 


VTIS  CRAfcl 
PTIC  TAB 

UnanncuocGd 
Jmstlf Ication. 


By - - 

Ion/ 


Availability  Codes 


OlSt 
.1 


111 


I'' 


\til  eo^or 
Spooial 


EXECUTIVE  SUMMARY 


The  problem  of  estimating  the  Doppler  frequency  of  a  received  RF  signal  is 
basically  the  problem  of  estimating  the  shift  in  frequency  of  the  signal  in  the  presence  of 
noise  generated  in  the  propagation  medium  and  in  the  measuring  systems.  This  paper 
establishes  bounds  for  the  best  possible  frequency  estimate  that  can  be  achieved  from  a 
finite  number  of  I  and  Q  samples.  The  latter  axe  obtained  by  demodulating  the  received 
signal  in  its  in-phase  (I)  and  quadrature  (Q)  components  which  are  then  discrete-time 
sampled. 

The  equations  describing  the  Cramer-Rao  (C-R)  lower  bounds  for  the  minimum 
variance  of  the  unbiased  estimates  of  the  frequency,  amplitude,  and  phaise  of  the  signal 
2ure  derived.  It  is  found  that  the  C-R  bound  for  the  frequency  estimate  is  inversely 
proportional  to  the  signai-to-noise  ratio  SNR,  the  square  of  the  sampling  time-interval 
T„  and  the  cube  of  the  number  of  samples  N  obtained.  The  minimum  variances  for  the 
coefficients  of  a  general  polynomial  phase  signal,  i.e,  chirp  signals,  are  also  derived.  The 
maximum-likelihood  (ML)  estimation  of  these  parameters  are  then  considered  since  the 
ML  estimator  is  the  minimum  variance  estimator.  The  analysis  of  the  ML  frequency 
estimator,  however,  reveals  that  at  low  SNR,  the  mean-square  error  (MSE)  of  the  ML 
estimate  increases  much  more  rapidly  than  the  expected  C-R  bound  due  to  the  higher 
probability  of  the  occurrence  of  an  outlier.  The  SNR  dependent  deviation  of  the  MSE 
from  the  C-R  bound  is  the  threshold  effect.  The  relationship  between  the  ML  estimation 
and  the  discrete  Fourier  transform  (DFT)  of  the  sampled  complex  I/Q  components  is  also 
shown. 


An  algorithm  for  the  implementation  of  the  ML  frequency  estimator  is  then 
considered.  It  consists  of  two  parts;  a  coarse  search  and  a  fine  search.  The  coarse  and 
the  fine  search  are  performed  using  the  fast  Fourier  transform  (FFT)  £«id  a  numerical 
procedure  such  as  the  bisection  method  or  the  secamt  method,  respectively.  Errors  due 
to  the  direct  estimation  of  the  frequency  parameter  by  the  FFT  algorithm  are  also  anad- 
ysed.  For  moderate  to  high  SNR  (above  the  threshold  effect),  the  dominant  error  in 
the  frequency  estimation  is  found  to  be  the  approximation  of  the  ML  estimate  by  the 
maximum  magnitude  DFT  point.  This  approximation  can  be  made  more  accurate  by 
increasing  the  number  of  DFT  points,  but  once  this  number  is  chosen,  the  MSE  of  the 
DFT  approximation  is  constant  and  is  independent  of  the  input  SNR. 

Some  numerical  examples  about  the  best  possible  frequency  estimates  that  can 


v 


be  achieved  for  coherent  and  noncoherent  pulse  radars  are  finally  given.  The  typical 
values  chosen  for  the  pulse  duration,  sampling  rate,  and  SNR  are  1  ns,  10®  Hz,  and  20  dB, 
respectively.  It  is  seen  that  in  the  case  of  coherent  pulses,  for  the  frequency  estimate  to 
have  an  accuracy  of  10  Hz  or  less,  the  number  of  pulses  required  is  around  50.  On  the 
other  hand,  in  the  noncoherent  case,  for  the  same  accuracy,  the  required  number  of  pulses 
is  in  the  order  of  a  100,000. 


VI 


TABLE  OF  CONTENTS 


ABSTRACT/RESUME .  iii 

EXECUTIVE  SUMMARY  .  v 

TABLE  OF  CONTENTS .  vii 

LIST  OF  FIGURES .  ix 

LIST  OF  TABLES .  xi 

1.0  INTRODUCTION  . • .  1 

2.0  COHERENT  RECEIVER .  6 

3.0  VARIANCE  BOUNDS  ON  PARAMETER  ESTIMATION .  7 

3.1  Cramer-Rao  Bounds  for  Frequency  and  Phase  Estimates .  8 

3.2  Cramer-Rao  Bounds  for  Polynomial  Phase  Modeled  Estimates  .  10 

4.0  MAXIMUM-LIKELIHOOD  ESTIMATION  .  14 

4.1  ML  Estimation  of  the  Polynomial  Coefficients .  15 

4.2  Relationship  with  the  Discrete  Fourier  Transform .  16 

4.3  Threshold  Effects  Due  to  the  Occurrence  of  Outliers .  20 

5.0  DETERMINATION  OF  THE  ML  FREQUENCY  ESTIMATE .  24  • 

6.0  SUMMARY  AND  CONCLUSION .  27 

6.1  Pulsed  Signals  .  28 

REFERENCES . REF-1 


vu 


LIST  OF  FIGURES 


Figure  1 :  Change  in  trajectory  of  the  target  emitter 

Figure  2:  Plot  of  the  range  R  versus  the  response  time  tr  for  a  trajectory  devi¬ 
ation  of  ^  =  5"  and  an  angular  resolution  of  0.3° . 

Figure  3:  Block  diagram  of  the  in-phase  and  quadrature  (I/Q)  demodulator  . 

» 

Figure  4:  Spectrum  of  1A(/)|  with  N  =  64,  fd  =  500  Hz,  /,  =  10  KHz,  and 

SNR  =  0dB . 

Figure  5;  Spectrum  of  \A{f)\  with  N  =  32,  fd  =  500  Hz,  /,  =  10  KHz,  and 
SNR  =  20  dB . 

Figure  6:  Probability  of  outlier,  q,  versus  SNR  for  various  vadues  of  iV . 

Figure  7:  Plot  of  RMS  error  of  the  ML  frequency  estimates  against  SNR,  with 

f,  =  10  KHz . 

Figure  8:  Plot  of  RMS  error  of  the  ML  frequency  estimation  against  SNR,  with 

U^lOOKHz . . 

Figure  9:  Overall  RMS  frequency  error  as  a  function  of  SNR  and  N  with  /,  = 

10  KHz  and  m  =  8 . 

Figure  10:  (a)  pulse  train  for  signal  frequency  /  >  (b)  pulse  train  for  signed 

frequency  f  <  \  and  /  <  PBF’,  (c)  pulse  train  for  signal  frequency 
fd<^^dfd>  PRF . 


LIST  OF  TABLES 


Table  1: 


Page 


Comparison  between  an  ESM  Doppler  Receiver  and  a  Monopulse  Radar 
System .  5 


% 


xi 


1.0  INTRODUCTION 


The  measurement  of  the  doppler  frequency  is  basically  the  measurement  of  the 
shift  in  frequency  of  the  received  signal.  In  the  case  of  a  radar  signal,  the  doppler  frequency 
shift  is  related  to  the  radial  component  of  the  relative  velocity  of  the  target  by  [9] 


/d  = 


2kf 

c 


(1) 


where  R  is  the  range  rate,  A  =  //c  is  the  transmitted  wavelength  of  the  radar  signal  and 
fi  is  the  doppler  shift  which  is  positive  for  decreasing  range  R.  For  a  passive  intercept 
receiver  such  as  the  Electronic  Intelligence  (ELINT)  or  the  Electronic  Support  Measure 
(ESM)  system,  the  doppler  shift  is  given  by  [12] 


f  -  ^ 


(2) 


where  R  is  the  relative  radisd  velocity  between  the  emitter  and  the  receiver. 


There  are  many  reasons  why  one  would  use  a  passive  system  such  as  the  ESM 
receiver  instead  of  an  active  radar  to  measure  the  doppler  frequency.  Among  them,  one 
obvious  reason  is  the  unavailability  of  a  radar  at  the  time  due  for  example  to  its  use  in 
tracking  other  targets.  More  importantly,  however,  the  ESM  system  has  the  advzmtage 
that  it  can  detect  an  emitting  target  much  earlier  than  a.  radar  because  of  the  one  way  wave 
propagation  instead  of  two.  Thus,  information  about  a  change  in  the  target’s  movement  or 
trajectory  can  be  made  available  much  sooner  using  an  ESM  system.  Since  the  intercept 
(ESM)  receiver  has  no  knowledge  of  the  transmitted  frequency  and  cam  only  measure 
the  received  signal,  it  cannot  calculate  the  shift  in  frequency  due  to  the  relative  motion 
between  the  emitter  and  the  receiver.  It  can,  however,  determine  the  doppler  shift  due  to 
the  change  in  the  radial  velocity  of  the  emitter.  This  can  be  useful  in  detecting  the  radial 
acceleration  or  the  change  in  trajectory  of  an  active  emitter.  As  am  exaunple,  suppose 
that  the  velocity  of  a  target  emitter  (missile)  is  Vi  m/s  amd  is  travelling  directly  towairds 
the  receiver  (ship  A)  ais  shown  in  Figure  1.  At  some  point  atlong  the  route  it  chamges 
trajectory  amd  turns  (fi  degrees  away  from  the  ship.  The  change  in  radial  velocity  due  to 
the  change  of  <f>  degrees  from  the  original  trajectory  is  K  —  K  cos  </>  =  Vi{\  —  cos^)  m/a. 
Hence,  the  shift  in  frequency  of  the  signal  at  the  receiver  is 


Figure  1:  Change  in  trajectory  of  the  target  emitter. 
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For  ^  =  5®,  A  =  0.03  m  and  Vi  =  300  m/s,  the  doppler  frequency  sliift  fd  is  then  equal  to 
38.05  Hz.  In  terms  of  the  angle  change  in  trajectory,  (^,  we  can  rewrite  equation  (3)  as 


(f)  =  cos 


where  fd  is  the  shift  in  frequency  and  /a  >  0.  Assuming  that  we  are  able  to  estimate  /a 
down  to  say  around  10  Hz.  With  the  same  values  of  K  and  A,  it  means  that  we  will  be 
able  to  detect  a  change  in  trajectory  of  at  least  2.56'^. 

Now,  we  will  discuss  the  advantage  of  using  an  ESM  system  with  the  capability 
of  measuring  the  doppler  frequency  with  an  accuracy  of  about  10  Hz  or  less  compared  to  a 
monopulse  radar  system  that  has  an  angular  resolution  of  0.3  degrees.  Referring  to  Figure 
1,  consider  the  first  scenario  where  the  missile  heis  locked  onto  the  ship  A  and  the  latter  is 
trying  to  transfer  the  lock  from  itself  onto  the  nearby  chaff  cloud  through  the  use  of  some 
electronic  countermeasures.  The  chaff  was  released  by  the  ship  in  an  attempt  to  seduce 
the  missile  and  as  such,  it  is  not  very  fai'  from  the  ship.  In  this  sctmario,  it  is  assumed  that 
the  distance  (d)  between  the  chaff  cloud  and  ship  A  is  around  300  metres.  To  simplify 
the  example  the  ship  and  chaff  cloud  are  in  a  line  perpenticular  to  the  initial  missile 
trajectory.  It  is  shown  that  for  a  receiver  which  can  measure  a  shift  in  frequency  down  to 
10  Hz,  it  can  detect  a  change  in  trajectory  of  at  least  2.56  degrees  for  V,  =  300  m/s  and 
A  =  0.03m.  Hence,  this  implies  that  the  receiver  will  be  able  to  detect  any  deviation  of 
the  missile  onto  the  chaff  within  the  range  i?enn  =  d/  tan  </>  =  0.3/  taji(2.56°)  =  6.703  Km. 

For  a  monopulse  ra^iar  system  with  am  amgular  resolution  of  ^  degrees,  the 
maodmum  range  R  for  which  a  change  in  trajectory  of  (/>  degrees  cam  be  detected  in  a  time 
tr  seconds  or  less  is 

^  trViSm(j>  _  trViSin<i> 
tan/?  ~  ^ 

where  for  smadl  vaJues  of  tamy^  ~  /?.  As  can  be  seen  from  the  above  expression, 
the  ramge  R  is  proportionad  to  the  required  response  time  tr  amd  the  velocity  V,  of  the 
target,  but  inversely  proportionad  to  the  amgulau:  resolution.  A  plot  of  range  (R)  versus 
the  response  time  {tr)  for  different  velocities  (V<)  is  shown  in  Figure  2  with  the  trajectory 
change  oi  <f>  =  amd  am  angular  resolution  of  /0  =  0.3".  In  the  first  scenairio  as  described 
above,  the  deviation  angle  </>  must  satisfy  the  equation  tarn  ^  =  d/(i?  +  trK  cos  (f>).  Solving 
simultameously  this  equation  with  that  of  (5),  we  obtain  the  maximum  range  Rradar  = 
5.12  Km  for  V,  =  300 m/s,  U  =  23,  0  =  0.3",  and  d  =  300m. 
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Figure  2:  Plot  of  the  range  R  versus  the  response  time  U  for  a  trajectory  deviation  of 
^  =  5®  and  an  angular  resolution  of  0.3**. 


In  the  second  scenario,  the  chaff  is  used  against  the  acquisition  mode.  That 
is,  the  chaff  is  placed  in  advance  before  the  missile  acquires  its  target  and  as  such,  it  is 
much  further  away  from  the  ship  than  in  the  seduction  mode.  In  this  case  the  distance 
d,  between  the  chaff  and  the  ship  is  tadcen  to  be  approximately  1  Km  as  shown  in  Figure 
1.  Using  this  value  of  d  and  with  the  other  vauriables  having  the  sarnie  vadues  as  before, 
the  matximum  ramge  R^m  for  which  the  ESM  doppler  receiver  cam  detect  a  change  in 
trajectory  is  22.37  Km.  Note  that  R^^m  is  proportional  to  the  distance  d  between  the 
ship  amd  the  chaff.  For  the  monopulse  radar  with  the  sarnie  angulau:  resolution  of  0.3 
degrees  and  maximum  response  time  of  2  seconds,  by  solving  the  previous  equations,  the 
maximum  range  Rrmdar  obtained  is  10.39  Km. 

Similarly,  in  the  third  scenario  where  the  missile  deviates  to  track  ship  B  instead 
of  ship  A,  the  doppler  receiver  can  detect  this  chamge  well  beyond  the  horizon.  With  the 
distamce  d  between  ship  A  and  B  being  around  3  Km^  and  with  all  the  other  variables 
remauning  the  same,  we  get  R^  =  67.10  Km  and  i?radar  =  18.13  Km.  The  results  au-e 
summarized  in  Table  1  for  the  three  different  scenarios. 

Ideally,  if  we  were  able  to  determine  any  shift  in  frequency  at  the  receiver,  we 
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Table  1:  Comparison  between  an  ESM  Doppler  Receiver  and  a  Monopulse  Radar  Sys¬ 
tem. 


Scenarios 

ESM  Receiver  with  Doppler 
measurement  of  an  accuracy 
of  10  Hz 

Radar  System  with  angular 
resolution  of  0.3°  and  a  required 
maximum  response  time  of  2  s 

1.  d  =  300  m 

=  6.710  Am 

Rradar  =  5.118  Km 

2.  d  =  1.0  Km 

Ream  =  22.367  Km 

Rradar  =  10.388  Km 

3.  d  =  Z.QKm 

Ream  =  67.099  Am 

Rradar  =  18.130  Km 

would  be  able  to  calculate  any  trajectory  change  of  the  target.  However,  this  is  usually 
not  the  case  due  to  the  presence  of  noise  in  both  the  propagation  medium  and  in  the 
systems.  This  paper  addresses  the  problem  of  estimating  the  frequency  of  the  signal 
^  which  is  corrupted  by  additive  white  Gaussian  noise.  The  signal  is  assumed  to  be  of  the 

form  x{t)  —  a(t)  cos(a>,<t  -)-  6),  where  the  amplitude  a(<),  the  phase  $,  and  the  angular 
frequency  are  all  assumed  to  be  unknown.  In  addition,  the  estimation  is  performed  on 
a  finite  number  of  discrete-time  observations. 

The  minimum  variance  or  mean-squaxe  error  (MSE)  of  the  unbiased  frequency 
estimation  is  derived  using  the  Cramer-Rao  (C-R)  bound.  This  enables  us  to  have  an 
idea  about  the  best  possible  accuracy  that  can  be  achieved  as  a  function  of  the  number  of 
samples  and  input  signal-to-noise  ratio  (SNR).  To  further  generalize  this  concept,  we  have 
•  also  included  signals  that  are  frequency  modulated  such  as  the  chirp  signals.  The  lower 

variance  C-R  bounds  for  the  estimation  of  the  coefficients  of  an  arbitreury  polynomial 
^  phase  signal  of  the  form  x{t)  =  a{t)cos{bo  +  bit  -f  62*^  +  •••■+  are  derived.  The 

maximum-likelihood  (ML)  parameter  estimation  of  the  deniudulated  signal  is  considered 
in  this  paper  because  it  is  known  that  for  large  values  of  SNR,  the  ML  estimator  is 
unbiased,  and  is  the  minimum-variance  estimator  whose  meem-square  error  is  given  by 
the  C-R  bound.  As  the  SNR  decreases,  however,  the  increase  in  the  probability  of  the 
outliers  of  the  frequency  estimate  occuring,  will  cause  the  MSE  of  the  ML  estimator  to 
increase  much  more  rapidly  than  the  C-R  bound.  The  tendency  of  the  MSE  to  start 
deviating  from  the  C-R  bound  as  the  SNR  decreases  below  a  certain  value  is  known  as 
the  threshold  effect.  Most  nonlinear  estimators  suffer  from  threshold  effects.  It  will  be 
also  shown  that  the  ML  estimation  of  the  frequency  parameter  is  related  to  the  discrete 
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Figure  3:  Block  diagram  of  the  in>phase  and  quadrature  (I/Q)  demodulator. 

Fourier  transform  (DFT)  of  the  complex  envelope  of  the  sampled  signal.  The  overall  MSE 
in  using  the  DFT  as  an  approximation  to  the  ML  frequency  estimation  is  then  derived. 


2.0  COHERENT  RECEIVER 


We  will  first  assume  that  the  incoming  RF  signal  from  the  emitter  is  coherent. 
It  can  be  a  continuous  wave  (CW),  a  trsdn  of  coherent  pulses,  or  any  transmitted  sign2d 
whose  carrier  is  modulated  by  a  waveform  a{t).  For  now,  we  will  further  assume  that  the 
initial  carrier  frequency  Wc  has  already  been  obtained  and  is  known.  Suppose  that  there 
is  a  shift  in  frequency  of  the  signal  due  to  the  doppler  effect.  The  incoming  signal  x{t)  is 
then  given  by 

x(t)  =  a{t)  cos[(u;c  +  +  ^o]  +  n{t)  (6) 

where  wj  =  2irfd  is  the  doppler  shift,  and  n{t)  is  the  additive  white  Gaussian  noise  with 
zero  mean  and  variance  <t^.  The  received  signal  x(t)  is  then  applied  to  the  in-phase  and 
quadrature  (I/Q)  demodulator  as  shown  in  Figure  3.  Note  that  in  practice,  the  RF  signal 
is  usually  down-converted  to  an  intermediate  frequency  (IF)  signal  before  it  is  applied  to 
the  I/Q  demodulator.  This  is  not  shown  here  for  brevity. 
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The  noise  n{t)  can  be  represented  as  a  band  noise  signal  in  terms  of  its  in-phase 
and  quadrature  components  ndt)  and  respectively,  where  ndt)  and  n,(<)  are  also 

Gaussian  with  zero  mean  and  variance  [4],  that  is, 

n{t)  =  nc(<)cos(u;ct)  -  n,(<)  sin(ci;cO-  (7) 

The  demodulated  in-phase  and  quadrature  baseband  signals  I{t)  and  Q{t)  are  then  given 
by 


/(<)  =  a(<)cos(a;d< -f  ^o)  +  nc(0  (8) 

Q{t)  =  a(<)sin(u;(it -f  0o)  +  (9) 

The  baseband  signals  I{t)  and  Q{t)  are  sampled  at  a  rate  of  --  jr  per  second.  The 
sampled  signals  obtained  axe  then  equal  to 

Ik  =  akCO3{udtk  +  0o)  +  nc  (10) 

Qk  =  aksm{udtk-\-9o)-\-n^,  (11) 


where  =  <o  +  kTs,  for  A:  =  0,  1,  2,  . . .  and  /<.,  Qk,  Ok,  Ug  and  n,  are  defined  as  I{tk), 
Q(^fc)>  <*(U)>  ndtk)  and  n,(tfc),  respectively.  The  noise  samples  tie  and  n,  are  statistically 
independent  Gaussian  random  variables  with  mean  zero  and  variance  tr*.  For  CW  radar 
signals  the  amplitude  a(t)  is  usually  a  constant  and  for  pulsed  r2uiar  signals  the  value  of 
a{t)  is  either  a  constant  or  zero,  neglecting  the  rise  and  fall  time  of  the  pulses.  In  the 
latter  case,  we  can  always  sample  the  signals  I{t)  and  Q{t)  such  that  the  value  of  Ck  is 
the  constant  but  not  equal  to  zero.  Thus,  we  can  replace  a*  simply  by  the  p2urameter  a^. 

3.0  VARIANCE  BOUNDS  ON  PARAMETER  ESTIMATION 

Before  the  estimations  of  the  parameters  such  as  the  frequency  u^d  and  the  phase 
Oo  are  made,  it  is  usually  very  useful  to  have  numbers  or  values  that  will  indicate  the 
best  possible  estimates  that  can  be  achieved  with  the  available  data.  Suppose  that  we 
have  N  sample  values  of  Ik  and  Qk,  i.e.,  A:  =  0,  1,  •  •  • ,  iV  —  1.  From  these  values  of  Ik 
and  Qk,  we  can  form  the  N  sampled  complex  envelope  Zk  =  Ik  +  jQk,  where  j  = 

Since  the  noise  samples  are  independent  Gaussian  random  variables  with  mean  zero  and 
variance  <t^,  the  joint  probability  density  function  (PDF)  of  the  random  variables  Z*, 
A;  =  0,  1,  •  •  • ,  iV  —  1,  or  of  the  random  vector  Z  =  [Zo  Zi  Z2  •  •  •  Zn-\],  is  the  product 
of  the  two  IV-variate  Gaussian  distribution  of  the  random  vectors  I  =[Iq  I\  I2  •  •  •  /iv-i] 
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and  Q  =  [Qq  Q\  Q2  •"  Qn-i]-  The  PDF  of  vector  Z  is,  therefore,  given  by  [7] 


p(Z;  w, a, 6)  =  +  iQk  - 

where  the  means  fit  and  i/k  are  as  follows: 


fik  =  acos(u;<fc  +  fl) 
Vk  =  asin(u;<A:  +  d). 


(12) 


(13) 

(14) 


As  can  be  seen  from  the  above  equation  (12),  the  PDF  of  the  random  vector  Z  is  a 
function  of  the  unknown  parameters  a;,  a,  and  0. 


3.1  Cramdr-Rao  Bounds  for  Frequency  and  Phase  Estimates 


Due  to  the  fact  that  the  received  signal  is  corrupted  by  noise  and  is,  therefore, 
random  in  nature,  we  can  only  obtain  estimations  of  the  signal  parameters.  Furthermore, 
since  the  values  of  the  estimates  are  dependent  on  the  observat’*'  data,  they  axe  also 
random  variables  and  as  such  have  a  mean  and  a  variance.  An  unbiased  estimator  is  one 
whose  mean  value  is  equal  to  the  true  value  of  the  parameter  and  a  minimum-variance 
estimator  is  one  whose  variance  is  less  than  or  equal  to  the  variance  of  finy  other  estimator. 
The  minimum- variance  vmbiased  estimator  is  unique,  but  it  may  not  always  exist.  Under 
certain  conditions,  it  can  be  shown  that  the  variance  of  any  unbiased  estimator  cannot 
be  less  than  a  particular  lower  bound  known  as  the  Cramer- Rao  (C-R)  boimd  [11,  10]. 

FVom  [11],  the  Cramer- Rao  boimd  for  an  unbiased  estimator  of  a  single  unknown 
parameter  a  is  derived  and  is  shown  to  be 


where  a  is  the  estimate  of  the  true  value  of  a.  For  the  variance  of  the  simultaneous 
estimates  of  two  or  more  unknown  parameters  ai,  02,  •  •  orn,  the  unbiased  C-R  bounds 
are  obtained  from  the  inverse  of  the  (n  x  n)  Fisher  information  matrix  F  whose  elements 
are  given  by 

^fdlnp(Z;  ai,a2,-*-,a„)  5  lnp(Z;  ai,a2,”-,a„)'( 

=  •®[ - §5; - - / 
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f 


f 
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_  lnp(Z;  Q!i,a2,--  ,Qn)l 

“  \  daidoj  J  ■ 


(16) 


If  we  denote  the  inverse  of  the  matrix  F  by  i.e.,  'I'  =  then  the  C-R  bound  for  the 
variance  of  the  individual  unbitised  estimate  is  obtained  from  the  diagonal  elements  of 
That  is, 

Var(Q!i)  =  <r|,  >  V’u  (17) 

where  i//u  is  the  tth  diagonal  element  of  the  matrix 

Now,  returning  to  the  case  where  the  PDF  of  the  random  variable  vector  Z  is 
as  given  in  (12)  and  the  unknown  parameters  oi,  oa,  and  03  are  w,  a,  and  9,  resp>  V, 
we  have  a  (3  x  3)  matrix  P  whose  elements  are 


v^«i/  V^«.7 


(18) 


where  fXk  and  i/k  are  as  shown  in  (13)  and  (14),  respectively.  Solving  for  the  previous 
expression  (18),  and  substituting  tt  =  to  +  ^Ti,  the  symmetric  matrix  P  is  reduced  to 


a^{Ntl  +  N{N-l)toTs  +  lNiN-l){2N-l)T^)  0  {Nt^  +  \N{N  -  1)T,) 

0  A  O 

{Nto  +  \N{N  -  1)Z)  0  a^N 


(19) 

Taking  the  inverse  of  P,  the  Cramer-Rao  bounds  for  the  unknown  parameters  u;,  a  and  9 
are  as  follows: 


Var(u;)  <7?  >  _  1) 

Var(a)  =  <r|  >  ^ 


(20) 

(21) 


Var(^')  =  (r|  > 


2(7^  {64  +  6(Af  -  l)toT,  +  iN-  1){2N  -  1)T^) 

a^T^N{N^  -  1) 


(22) 


Note  that  unlike  the  frequency  u  and  the  amplitude  a,  the  variance  bound  on  the  phase 
9  is  dependent  on  the  time  at  which  the  first  sample  is  taken,  i.e.,  to-  These  results 
were  similarly  obtained  by  Rife  [8].  By  differentiating  equation  (22)  with  respect  to  to 
and  setting  the  result  to  zero,  it  is  found  that  the  minimum  bound  is  achieved  when 
to  =  — |(Af  —  1)T,  and  is  equal  to  <T^/{a^N).  In  terms  of  the  signal-to-noise  ratio  which 
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is  defined  as  SNR  =  a^/(2a^),  the  Cramer-Rao  bound  for  the  variance  of  the  frequency 
estimate  is  <t?  >  6fj/[{SNR)N{N'^  —  1)],  where  /,  =  l/T,  is  the  sampling  rate. 


3.2  Cramer-Rao  Bounds  for  Polynomial  Phase  Modeled  Estimates 

We  can  further  generalize  the  incoming  signal  x{t)  as  an  arbitrary  polynomial 
chirp  signal  given  by 

x(<)  =  a  cos  (bo  +  (oJc  +  6i)t  +  62^^  +  •  •  •  -I-  4-  n{t)  (23) 

where  a,  bo,  bi  bp  are  the  unknown  parameters  to  be  estimated.  That  is,  there  is  now 
a  modulation  on  the  frequency  or  phase  of  the  signal.  Note  that  for  p  =  1,  the  signal 
is  basically  the  same  as  that  in  the  previous  case  where  the  sinusoidal  has  a  constant 
frequency  equal  to  61.  For  p  =  2,  we  have  a  linear  chirp  or  linear  frequency  modulation 
(FM)  signail,  i.e.,  its  frequency  varies  linearly  as  a  function  of  time  and  the  slope  is  given 
by  63.  Higher  order  chirp  signads  can  be  obtained  with  p  >  3.  As  before  using  the  same 
coherent  receiver  as  that  described  in  section  2,  we  have 

/(<)  =  a  cos(bo  +  bit  +  bjt^  +  •  •  •  +  V**)  +  (24) 

Q(t)  =  a  sin(bo  + bit  +  b^t^ +•••+ bpt'")  +  nf(t).  (25) 

The  above  in-phase  and  quadrature  signads  aue  sampled  into  the  discrete  samples  Ije, 
and  Qk,  k  =  0,  1,  . . .,  iV  —  1,  respectively.  Defining  the  complex  samples  Zk  =  h  + 
jQkt  the  joint  probability  density  function  of  vector  Z,  given  the  parameter  vector  5  = 
[a  bo  bi  ...  6p],  i.e.,  p(Z;  a),  is  the  same  as  that  expressed  in  (12),  except  that  the  meams 
/ik  and  Vk  axe  now  given  by 

fik  =  a  cos(6t)  4-  bitk  4-  hjtfc  4-  •  •  •  4-  bptl)  (26) 

i/k  =  a  sin(6b  4-  bitk  4-  62^*  4*  •  •  •  4-  bptl)  (27) 

It  is  to  be  noted  that  fik  amd  Vk  are  now  functions  of  the  (p  4-  2)  unknown  parameters 
a,  bo,  hi, . . . ,  bp.  Similar  to  the  previous  example,  the  Cramer-Rao  bounds  for  the  vauiotis 
parameters  can  be  obtained  from  the  inverse  of  the  ((p  4-  2)  x  (p  -f-  2))  Fisher  informa¬ 
tion  matrix  F  whose  elements  7,^  are  given  by  equation  (18),  and  where  the  unknown 
parameters  a<,  t  =  1,  2,  . . . ,  p  4-  2,  are  a,  60,  61,  ... ,  bp,  respectively.  Solving  for  7i,, 
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t,  j  =:  1,  2,  •  •  • ,  p  +  2,  we  have 


t 

where  the  elements  Si  are 


(  ^ 

0 

0 

0 

0 

5-0 

••• 

Sp 

a* 

0 

5, 

52  ••• 

Sp+i 

.  0 

^P 

•^p+i 

S2p 

N-l 

N-l 

Si  = 

k=0 

(tkY  = 

E(<o  +  A:7;)‘. 
fc=o 

(28) 


(29) 


The  matrix  F  can  now  be  inverted  and  the  diagonal  elements  of  F”^  yield  the  C-R  lower 
bounds  for  the  variance  of  the  unbiased  estimates  of  each  of  the  parameters  a,  bo,  bi,  '  •  • ,  bp. 
Note  that  the  inverse  of  the  matrix  F  can  be  computed  by  reducing  it  to  the  (p+1)  x  (p+1) 
matrix  P,  where 

Si  Si  •  •  •  Sp 

Si  Sz  • '  • 

Si  S3  S4  •  •  •  Sp+2 


(  So 
Si 


F'  =  ^ 


Sp  5, 


p+l  •^p+2 


(30) 


; 


Once  the  inverse  of  F'  is  obtained,  then  F  *  is  simply  the  augmented  matrix  of  (F')~^  as 
shown  below. 

0  \ 


F-'  = 


/  0 
0 


i  [  (ro-  ] 


k  0 


(31) 


What  is  required  now  is  to  calculate  (F')~^.  The  matrix  F'  has  a  Hankel  type  structure 
and  hence,  can  be  efficiently  computed  using  such  recursive  procedure  as  the  Levinson- 
Durbin  algorithm  [3].  We  will  now  consider  a  (4  x  4)  matrix  F'  which  implies  a  third 
order  polynomial  phase  modeled  of  the  signal,  i.e.,  x(t)  =  a  cos(6o  +  bit  +  bit^  +  63<^).  The 
inverse  matrix 


•  'fill  'fu  V’n  <^14  ^ 

0^  V’ai  <P22  ^*23  <^24 

D  ^32  (p33 

<P*l  V^42  V?43  1^44  J 


(32) 
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where  the  determinant  D  is  given  by 


D:=  St-  3525^54  +  SlSl  +  2SiS:iSl  -  SqS^  +  2SIS3SS  -  2SiSlSs  -  2S,S2S4Ss 
+25o535455  +  SlSl  -  S0S2SI  -  S^Ss  +  23^823386  -  S^SlSa  -  SfS^Se 

+S0S2S4S6  (33) 

and  the  elements  (pi  ,j,  i,  j  =  1,  2,  3,  4,  are 


<P\i  = 

<Pl2  —  *P21  = 
<Pl3  —  fPSl  = 
9l4  =  <^41  = 
(P22  = 

<^23  =  <P32  = 
<^24  =  9>42  = 
<P33  = 

<P34  —  <P43  = 
¥>44  = 


823486  -  53  +  283S4SS  -  S2SI  -  SjSs 

S3SI  -  SiSs  -  82848s  +  8i8i  +  S'25356  -  5t5456 

525i  -  5|54  +  52535s  -  818483  -  8l8e  +  S1S3S6 
SI  -  2S2S3S4  +  SiSl  +  SlSs  -  818383 
808486  -  SIS4  +  2S2S3S3  -  SoSl  -  8186 

83  —  S2S3S4  —  5i5355  +  5o5455  +  5i5256  —  S0S3S6 

SIS4  -  5251  +  S1S3S4  -  5o5|  -  5i5255  +  S0S3S3 

S0S2S6  -  S2SI  +  25i5354  -  5o5i  -  5?56 

8^83  -  SiSt  -  S1S2S4  +  S0S3S4  +  SIS3  -  S0S2S3 

S0S2S4  -  5|  +  25i5253  -  5o5|  -  5^54 


For  simplicity  it  is  assumed  that  the  first  sample  is  taken  at  to  =  0,  which  means  that 
tif  =  kTt-  In  this  case,  the  elements  Si  of  the  Fisher  matrix  F  becomes  5,-  = 

'n  For  the  (4  x  4)  matrix  F',  the  elements  So  to  So  are  as  follows:(refer  to  [1]) 

So  =  N 

5.  =  Mr. 


s,  = 


N(N  -  l)(2jV  -  1) 


iV’(JV  -  D’  ^3 
S3  -  - j - T, 

„  fUN  -  l)(2N  -  1)(3JV’  -  3rtr  -  1)  ^ 
=  - 30 - 


S3  = 


N^N  -  l)^(2iV^  -  2jV  -  1)  ^5 


^  N{N-l){2N-l){3N*-&N^  +  W  +  l)^ 

=  - 42 - ^ 

Note  that  we  are  only  interested  in  the  diagonal  elements  of  (F')“*.  The  Cramer- Rao 
lower  bounds  on  the  parameter  estimates  d,  60,  61,  62,  ^  (the  superscript  '  denotes  the 
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estimator  of  the  given  parameter)  are  then  given  by 


Var(a)  =  >  — 


Var(M  =  <  > 


Var(M  =  .?  >  = 

Var(^)  =  <  > 


4(2Af  -  l)(iV2  -  iV  +  3) 
{SNR)NiN  +  1)(A^  +  2){N  +  3) 

100(6A^‘‘  -  27 +  42N^  -  30N  +  11) 
{SNR)  N{N^  -  l)(iV2  -  4){N^  -  9)  Tj 

180(9^  -  13)(2^  -  1) 

{SNR)  N{N^  -  l)(iV2  -  4)(JV2  -  9)  Tf 

1400 

{SNR)  N{N^  -  1){N^  -  4){N^  -  9)  Tf 


(36) 

(37) 

(38) 

(39) 

(40) 


where  we  have  defined  the  signal-to-noise  ratio  SNR  =  a? (20^ .  These  expressions  can  be 
simplified  if  we  assume  that  the  number  of  samples  N  is  large,  i.e.,  >■  1,  In  this  case, 

the  value  of  Si  reduces  to 


Si 


N-\ 

=  r.  E 


ksO 


(41) 


Substituting  the  above  equation  in  (33)  and  (34),  we  have  the  approximate  values  of  the 
variance  lower  bounds  of  the  parameters’  estimates  6o,  to  63  as 


Var(M  =  <T^  > 


8 


{SNR)N 
600 


Var(6i)  <Tg^  >  (5^^) 

V!a{b2)  =  <Tl  > 

^  {SNmrnr 


(42) 

(43) 


(44) 

(45) 


FVom  the  expression  of  the  lower  variance  bounds  for  the  estimation  of  the  polynomial 
phase  coefficients,  it  is  seen  that  all  the  variances  are  inversely  proportional  to  the  signal- 
to-noise  ratio  as  expected.  In  addition,  it  can  be  observed  that  the  variance  bound  of  the 
polynomial  coefficient  is  inversely  proportional  to  twice  its  polynomial  power  of  the  time 
duration  r  =  NT,  of  the  signal,  where  N  is  the  number  of  samples  and  T,  is  the  sampling 
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interval.  For  example,  the  estimation  of  the  coefficient  62,  the  second  order  polynomial,  has 
a  variance  which  is  inversely  proportional  to  t*  =  (NT^y.  This  means  that  the  higher  the 
order  of  the  polynomial  coefficient,  the  more  number  of  samples  or  a  longer  time  duration 
of  the  signal  is  required  to  achieve  a  certain  accuracy  in  the  estimation  of  that  coefficient. 
It  is  to  be  noted  that  the  exact  expressions  in  (37)-(39)  depend  on  the  initial  sampling 
time  <0-  Their  minimum  bounds  are  usually  achieved  when  to  =  T,. 

4.0  MAXIMUM-LIKELfflOOD  ESTIMATION 

• 

Having  determined  the  lower  bounds  for  the  mean-square  errors  of  the  parame¬ 
ters’  estimates,  we  will  now  try  to  obtain  these  estimates  using  the  maximum- likelihood 
(ML)  method.  The  ML  estimator  maximizes  the  probability  density  function  p(Z;  a) 
with  respect  to  the  vector  a,  where  the  elements  of  a  are  the  parameters  to  be  estimated 
and  Z  is  the  observed  sample  data.  It  is  to  be  noted  that  the  maximum  of  p(Z;  a)  occurs 
at  the  same  point  as  the  maximum  of  lnp{Z]  a).  Hence,  the  maximization  of  the  PDF 
p(Z]  a)  is  equivalent  to  the  maximization  of  lnp(Z;  a).  We  will  first  consider  the  case 
of  the  simple  sinusoid  waveform.  From  (12),  we  can,  therefore,  maximize  the  expression 

lnp(Zi  =  C  -  5^^’  [(A  -  pO'  +  (<?*  -  ■'.)’]  (46) 

where  C  is  a  constant,  independent  of  the  parameters  u;,  a,  and  0.  By  noting  that 
the  summations  of  and  Ql  over  k  are  constcints  once  an  observation  has  been  made, 
and  eliminating  these  and  some  other  constants,  we  obtain  a  simplified  expression  to  be 
maximized  as 

N-l 

L  =  2Y,ihfik  +  Qkn)-Na^  (47) 

k=0 

where  pk  and  i/j^  are  as  defined  in  (13)  euad  (14),  respectively,  and  +  ^1)  = 

a*  =  To  find  the  maximum  of  L  we  first  observe  that  we  can  rewrite  (/*p*  -h 
QkVk)  in  terms  of  the  real  part  of  the  complex  envelope  Substituting  tk  by 

<0  +  kT,  and  if  normalizing  the  result  by  the  number  of  points,  N,  we  get 


Lo  =  2a»{>l(u;)e--'^e--''*'‘»}  - 

(48) 

with 

k=0 

(49) 
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where  3ii[*]  denotes  the  real  part  of  [•]. 

Since  the  p2Lrameters  w,  a  >  0,  and  9  are  all  unknown,  it  can  be  seen  that  Lo  is 
maximized  over  5  if  =  |>l(w)l.  This  occurs  when  ^  =  arg{A(a;)e~^‘^‘®}, 

where  arg[-]  is  defined  as  the  argument  or  phase  of  [•].  That  is, 

maxio  =  2a|A(u;)|  —  a^.  (50) 

6 

K  u^niM  is  the  value  of  u  which  maximizes  |A(u;)|,  then 

maxlo  =  2a|A(u;m»x)|  -  (51) 

9,0/ 

In  order  to  find  the  value  of  a  which  maximizes  Lo,  the  above  equation  can  be  differentiated 
with  respect  to  a  and  the  result  set  to  zero.  Performing  this,  we  obtain  a  =  |A(wn»ax)|> 
and 

max  To  =  |A(u;max)|^-  (52) 

9,0/ ,(1 

Hence,  the  m2ocimum-likelihood  estimates  for  Wd,  a^,  and  6o  denoted  by  tDj,  do,  and  9o, 
respectively,  are 


u)d  =  max|A(u;)|  (53) 

So  =  \A{C}d)\  '  (54) 

9o  =  arg{A(wd)e-^"«'‘»}.  (55) 

It  is  to  be  noted  that  the  optimization  of  the  function  Lo  C2ui  be  performed  as  described 
above,  because  Lo  is  separable  in  terms  of  its  variables  u,  a,  and  9. 


4.1  ML  Estimation  of  the  Polynomial  Coefficients 

For  the  general  case  of  a  sinusoid  having  a  polynomial  phase  model  as  described 
in  (23),  the  maximization  of  lnp(Z;  a),  where  the  parameter  vector  a  =  [a  6o  6i  •  •  •  6^], 
can  be  similarly  reduced  to  the  expression  L  in  (47),  but  with  /ifc  and  Vk  given  by  (26) 
and  (27),  respectively.  In  this  case,  (/*/<*  +  Qk^k)  can  be  rewritten  in  terms  of  the  real 

part  of  the  complex  envelope  Z*  aexp{-j(6o  +  hitk  +  htl  ^ - h  bptl)}.  As  before,  if  we 

assume  to  =  0,  then  tk  =  kTa,  and  after  normalization  by  N  (the  number  of  points),  the 
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expression  to  be  maximized  is 


Lo  =  2a  5(6,,  62,  •  •  • ,  6p)}  - 


(56) 


where  the  function 

2  N~l 

5(6i,  62,  •  •  • ,  6p)  =  —  ^Zk  exp  {-j[bi{kT,)  +  +  >  •  •  +  bp{kT,Y])  .  (57) 

*=0  ' 

Similar  to  the  previous  case,  Lo  is  maximized  over  bo  when  bo  =  ajg{B(6,,  62,  •  •  • ,  6^)}. 
With  this  value  of  60,  we  have 


m^Io  =  2al5(6i,  62,  •••,  6p)| -a^  (58) 

K  •  j  6p  are  the  values  of  6,,  62,  •  •  • ,  6p  that  maximizes  15(6,,  62,  •  •  • ,  6p)|  and 

which,  therefore,  also  maximizes  Lo^  then 


max  Lo  =  2al5(6i,  62, 

»  “‘lOp 


6,)i  - 


(59) 


The  maximization  of  the  above  expression  (59)  over  the  parameter  a  occurs  when  a  = 

A  A  A 

15(6i,  62,  •  •  • ,  5p)!  and  hence,  the  maximum  of  Lo  over  all  the  parameters  (a,  60, 6,,  •  •  • ,  bp) 
is 


max 

0.60,61  •••,6p 


Lo  —  |B(6i,  62, 


(60) 


As  can  be  seen,  this  result  is  a  simple  generalization  of  the  result  in  (48)-(52).  To 
solve  for  the  maximum  likelihood  estimation  of  the  polynomial  phase  signal,  requires  the 
maximization  of  the  function  |5(6i,  62,  •  •  • ,  6p)j  over  a  p-dimension2il  space,  which  is  a 
quite  difficult  problem  for  large  values  of  p.  From  here  on,  only  the  simple  case  where  the 
signal  is  a  sinusoid  having  a  fixed  frequency,  i.e.,  p  =  1,  will  be  considered. 


4.2  Relationship  with  the  Discrete  Fourier  Transform 


From  the  definition  of  the  discrete  Fourier  tr<insform  of  the  sampled  complex 
envelope  {Zk)  [5],  and  normalized  by  the  number  of  points  N,  we  have 


2  JV-i 


for  /  =  0,  1,  2,  ■  •  • ,  TV  -  1. 


(61) 
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Thus,  from  the  expression  of  A{u})  in  (49),  it  can  be  seen  that 

>1;  =  A  (^)  =  A  (5^)  tor  1  =  0,  1,2,. ...W-l  (62) 

where  a;,  =  27r/Tj  =  2irf,.  The  values  of  Ai,  therefore,  represent  the  sampling  points 
of  the  continuous  spectrum  >l(w),  with  the  sampling  intervals  equal  to  u;,/iV.  Note  that 
A(a;)  repeats  itself  at  every  w,  interval,  that  is,  24(u;)  is  periodic  with  period  Hence,  to 
find  the  maximum  of  |A(a;)|,  we  need  to  search  only  in  the  interval  between  0  and  u;,.  In 
addition,  if  we  have  an  a  priori  knowledge  about  the  range  or  the  maximum  value  of  the 
frequency  to  be  estimated  we  can  reduce  the  search  space  by  choosing  a  sampling 
rate  w,  close  to  the  largest  possible  value  of  wj. 

Figures  4  and  5  show  the  spectra  of  |A(/)|  (/  =  ^),  with  N  =  64,  fd  =  500  Hz, 
f,  =  \QKHz,  and  SNR  =  OdB  in  Fig.  5  and  iV  =  32,  fd  =  500  ifc,  /,  =  \QKHz, 
and  SNR  =  20  dB  in  Fig.  6.  The  crosses  along  the  curves  of  i  he  figures  represent  the 
magnitudes  of  the  N  points  DFT  of  A{f),  that  is,  the  {|A;i|  points,  /  =  0,  1,"-Ar  —  1. 
The  value  of  /  at  which  the  maximum  point  of  |24(/)|  occurs  is  the  ML  estimation  of  fd, 
i.e.,  fd,  and  the  value  of  \A(f)\  at  that  point  is  the  ML  estimation  of  Oo,  i.e.,  do.  In  the 
case  of  Figure  4,  the  values  of  fd  and  do  obtained  are  fd  =  507.23  ife  and  do  =  0.961, 
whereas  those  of  Figure  5,  fd  =  499.27  Hz  and  do  =  0.985.  From  the  N  DFT  points  in 
Figure  4,  the  maximum  is  IA4I  which  has  a  value  of  0.86  and  which  is  at  468.75  Hz  and 
in  Figure  5,  the  maximum  is  j/lal  with  a  vjdue  of  0.75  and  at  02b  Hz. 

Comparing  the  two  plots,  it  can  be  observed  that  having  a  larger  value  of  N 
causes  the  width  of  the  peaks  (lobes)  in  the  spectrum  to  become  smaller,  that  is,  the  width 
is  inversely  proportional  to  NT,.  As  the  signal-to-noise  ratio  increases,  the  sidelobes’  peaks 
decrease  and  the  values  of  fd  and  do  approach  the  exact  vadues  as  expected.  It  is  to  be 
noted  that  if  no  noise  is  present,  |A(u;)|  will  be  a  sine  function  which  is  symmetric  about 
the  doppler  frequency  and  which  has  a  period  of  u;«.  The  global  maximum  point  will 
then  occur  at  exactly  uid  and  will  have  an  exact  amplitude  of  a,,.  In  the  presence  of  noise, 
however,  the  sine  function  of  |A(a;)(  is  distorted  with  the  sidelobes’  peak  values  becoming 
larger  and  the  global  maximum  point  being  shifted  away  froin  the  true  value  of  wj.  The 
difference  between  the  maximum  DFT  point  to  the  true  maximum  of  |A(u;)|  depends  on 
M,  the  number  of  DFT  points  used  in  the  computation,  but  is  independent  of  the  SNR. 
The  larger  the  value  of  M,  the  more  accurate  is  the  DFT  approximation. 
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Figure  4;  Spectrum  of  |A(/)|  with  N  =  64,  /d  =  500  Hz,  /,  =  10  KHz,  and  SNR  -  0  dB. 
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Figure  5;  Spectrum  of  |A(/)|  with  N  =  32,  fd  =  500  /,  =  10  KHz,  and  SNR 


20  dB. 


4.3  Threshold  Effects  Due  to  the  Occurrence  of  Outliers 


In  this  section  we  will  consider  only  the  case  of  the  ML  estimation  of  the  fre¬ 
quency  parameter.  As  stated  previously,  for  large  SNR,  the  mean-square  error  of  the  ML 
estimation  can  be  obtained  from  the  Cramer-Rao  bound  given  in  (20).  However,  as  the 
SNR  decreases,  the  probability  of  occurrence  of  an  outlying  frequency  component  whose 
magnitude  is  the  largest,  increases.  An  outlier  is  defined  as  that  frequency  component 
which  is  not  in  the  main  lobe  containing  the  trcinsmitted  signal  frequency  but  which  has 
a  magnitude  larger  than  the  main  lobe.  When  the  maximum  peak  frequency  component 
(outlier)  occurs  far  removed  from  the  exact  frequency  mjj,  severe  errors  occur  in  the  ML 
frequency  estimate.  It  is  the  nature  of  these  large  errors  caused  by  the  presence  of  outliers 
that  will  be  examined  below. 

For  simplicity  of  the  analysis,  it  will  be  assumed  that  an  outlier  can  occur  any¬ 
where  between  zero  and  u;,.  That  is,  its  probability  density  function  is  approximately 
imiform.  In  this  case,  the  variance  or  the  MSE  due  to  an  outlier  is  U)l/12.  The  MSE  of 
the  ML  estimate  for  any  value  of  SNR  is,  therefore,  given  by 


MSE  =  P(outlier)  •  (MSE  due  to  outlier)  -f  P(no  outlier)  ■  (MSE  due  to  no  outlier) 

«  9  •  +  (1  -  9)  •  ^(5]v;2);/v(Ar2  -  1) 

where  q  is  the  probability  of  an  outlier  and  in  the  absence  of  any  outlier,  the  MSE  is 
taken  from  the  C-R  bound.  The  expression  for  q  was  derived  in  [8]  by  defining  the 
random  variable  Xi  —  |Aj|,  for  /  =  0,  1,  •  •  • ,  iV  —  1,  Xi  has  a  Rayleigh  distribution  when 
it  is  any  frequency  component  other  than  the  frequency  of  the  signal,  i.e.,  the  outlier,  and 
has  a  Rician  distribution  when  it  is  the  sign2d  frequency  component.  The  expression  of  q, 
thus  obtained,  was  shown  to  be 


q  =  1 


Nx 


Nx  iV(x^-l-a*)l  ,  ^Nax\  , 

75- “Pi - ^ 


(64) 


where  /o(')  is  the  modified  Bessel  function  of  the  first  kind  and  the  signal-to-noise  ratio 
(SNR)  =  a^/(2<r^)  is  as  defined  before.  For  large  values  of  N,  the  summation  representa¬ 
tion  in  (64)  is  not  easily  computed  because  of  the  factorial  term  which  can  get  extremely 
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Figure  6:  Probability  of  outlier,  g,  versus  SNR  for  various  values  of  N. 
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Figure  7:  Plot  of  RMS  error  of  the  ML  frequency  estimates  against  SNR,  with  /,  = 
IQ  KHz. 

large.  In  these  cases,  the  numerical  evaluation  of  the  integral  representation  is  more  prac¬ 
tical.  The  plots  of  q  versus  SNR  for  different  values  of  N  2ure  shown  in  Figure  6.  With 
these  values  of  9,  the  mean-square  error  of  the  maximum-likelihood  frequency  estimation 
given  by  (63)  can  be  calculated.  The  root  mean  square  (RMS)  error  of  the  ML  estimation 
is  then  obtained  by  taking  the  square-root  of  the  MSB.  The  RMS  frequency  errors  are 
plotted  against  the  input  SNR  for  different  values  of  N  in  Figures  7  and  8,  where  the 
sampling  rate  /,  j=  a;,/(2?r)  are  \QKHz  and  lOOiifflr,  respectively.  As  can  be  seen  from 
the  plots  of  Figxires  7  and  8,  depending  on  the  values  of  N  and  /„  the  threshold  effect 
takes  place  somewhere  below  QdB.  The  straight  lines  in  the  graphs  represent  approxi¬ 
mately  the  C-R  bounds.  For  a  given  value  of  N,  SNR,  and  /„  it  can  be  easily  determined 
if  the  ML  estimation  is  being  operated  above  or  below  the  threshold  effect  directly  from 
the  plots  in  Figures  7  and  8.  For  example,  if  iV  =  128,  SNR  =  20  dB,  and  /,  =  10  KHz, 
then  from  the  graph  of  Figure  7,  the  ML  estimation  is  above  the  threshold  and  the  RMS 
error  of  the  frequency  estimate  is  approximately  equal  to  the  C-R  bound  which  is  1.69 
rad/s  =  0.26  Hz. 
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Figure  8:  Plot  of  RMS  error  of  the  ML  frequency  estimation  agmnst  SNR,  with  /, 
imKHz. 


23 


5.0  DETERMBVATION  OF  THE  ML  FREQUENCY  ESTIMATE 


It  is  shown  in  Section  4,  that  the  ML  frequency  estimation  of  the  signal  is  the 
angular  frequency  u  that  maximizes  the  magnitude  of  A(u;),  where  A{u})  is  as  given  in 
(49).  To  find  m^|/l(c4j)|,  we  could  differentiate  the  function  \A{uj)\'^  =  A{u>)A*{u})  with 
respect  to  u)  and  set  the  result  to  zero.  This,  however,  does  not  usually  work,  since  |  A(u;)| 
or  |A(u;)|^  is  not  a  monotonic  function  of  u;  but  has  numerous  local  maxima  and  minima. 
In  order  to  determine  the  global  maximum  of  lA{u;)|,  we  first  require  a  rough  or  coarse 
estimate  of  where  it  is  located.  Once  this  is  obtained,  a  fine  search  can  then  be  performed 
to  find  the  exact  location  of  the  global  maximum  point. 

The  coarse  estimate  of  where  Wm.v  =  mjoc|A(u;)|,  can  be  made  directly 
from  the  computation  of  the  DFT  of  the  complex  data  {Zk},  k  =  0,  1,  •••  N  —  I  as  v/as 
done  in  [6].  With  the  value  of  lug/M  for  which  |A/1  is  maximum,  as  2m  initial  guess,  the 
bisection  or  the  secant  method  can  be  used  to  determine  Wmax  at  a  desired  accuracy.  To 
do  this,  we  need  to  evaluate  the  derivative  function 

=  =2»|A-(a.)^(A(u.)l|  (65) 

where  A*(w)  is  the  complex  conjugate  of  A(u;)  and  is  the  real  part  of  {•}.  Both 
the  bisection  and  the  secant  methods  use  an  iterative  procedure  to  compute  and  at 
the  maximum  point,  =  0.  The  initial  interval  in  which  the  numerical  algorithm 

is  applied  is  usually  taken  to  be  the  distance  between  two  consecutive  DFT  points,  i.e., 
u;,/M,  where  M  >  N.  This  will  ensure  that  there  exists  at  most  one  peak  frequency 
component  in  this  interval,  since  the  widths  of  the  peaks  £U‘e  inversely  proportional  to 
NT,. 

For  both  numerical  procedures,  we  need  to  evaluate  F(u;„+i)  at  each  iteration. 
In  the  case  of  the  bisection  method,  Wn+i  =  5(cn +dn),  where  c„  emd  d„  are  the  end-points 
of  the  interval  after  n  iterations  and  for  the  secant  method,  Un+i  is  calculated  from 

Wn+l  =  -  F{un)  — y  (66) 

The  secant  method  usually  has  a  faster  rate  of  convergence  than  the  bisection  method.  If 
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we  assume  that  after  i  number  of  iterations  the  numerical  procedure  is  stopped,  then  the 
computational  complexity  is  of  the  order  of  Z(N  complex  multiplications  and  2£N  complex 
additions  due  to  the  evaluation  of  F(u>„).  To  determine  the  initial  coarse  estimation,  we 
compute  an  M-point  DFT,  where  M  >  N.  This  can  be  performed  using  the  iast  Fourier 
transform  (FFT)  which  requires  A/logjM  complex  multiplications  and  additions  when 
A/  is  a  power  of  2.  For  moderate  to  high  SNR,  the  coarse  search  can  usually  be  done 
with  M  —  N.  The  overall  computational  complexity  using  this  method  of  coarse  and  fine 
search  is,  therefore,  of  the  order  of  A^(log2  N  +  £)  complex  multiplications  and  additions. 

An  alternative  to  the  previous  search  method  is  to  directly  find  an  approximate 
value  of  utxamx  by  having  a  very  fine  sampling  of  the  continuous  spectrum  of  |A(a;)|.  This 
can  be  achieved  by  computing  a  large  number  of  DFT  points.  That  is,  an  M-point  FFT 
is  performed  on  the  complex  data  {Zk},  k  =  0,  1,  •  •  • ,  /V  —  1,  where  A/  is  a  multiple  of  N, 
i.e.,  M  =  mN,  where  m  >  1.  In  this  case  the  computational  complexity  is  of  the  order 
of  Mlog^M  =  mN(logjN  +  logjm)  complex  multiplications  and  additions.  Note  that 
for  M  >  N,  the  N  points  data  {Zk}  is  padded  with  trailing  zeros  to  length  M  before  the 
DFT  or  FFT  is  performed.  The  approximate  value  of  u;m«v  is  then  given  by  •  w,/A/, 
where  =  mjw  |A/|  and  /  =  0,  1,  •  •  • ,  Af  —  1. 

The  exact  value  of  can  be  assumed  to  be  situated  anywhere  between 
and  its  neighboring  DFT  point.  This  implies  that  u;m«v  is  a  random  V2triable  that  is 
uniformly  distributed  in  the  interval  of  u;,/Af.  Hence,  the  MSE  associated  with  the 
approximation  of  ac,..  by  •  a;,/A/  is  equal  to  a;J/(12A/^).  The  overall  mean- square 
error  in  determining  the  signal  or  doppler  frequency  Ud  is,  therefore,  equal  to  the 
variance  of  the  ML  estimation  of  Ud  in  the  presence  of  noise  as  given  by  (63)  plus  the 
variance  of  approximating  UmMx  hy  the  maximum  magnitude  DFT  point.  That  is, 


2 


+  (1-9) 


+ 


u;: 


2ir^iSNR)N{N^  -  1)  12M^ 


(67) 


with  M  =  mN,  and  m  >  1.  For  moderate  to  high  signal- to-noise  ratio,  the  variance  of 
the  ML  estimation  can  be  approximated  by  the  C-R  bound  and  the  overall  MSE  can 
be  determined  from 

^  ^  2ir^SNR)N{N'^  -  1)  12m2iV2’ 

The  overall  root-mean-square  error  e  is  then  obtained  by  taking  the  square-root  of  e^. 
Depending  on  the  values  of  m  and  N  used,  it  is  seen  from  the  above  expression  (68)  that 
the  MSE  due  to  the  use  of  the  maximum  magnitude  DFT  point  as  an  approximation  to 
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Figure  9:  Overall  RMS  frequency  error  as  a  function  of  SNR  and  N  with  f,  =  10  KHz 
and  m  8. 

the  signal’s  frequency  (second  term)  can  be  much  larger  than  the  inherent  MSE  of  the  ML 

frequency  estimation  (first  term).  At  a  low  SNR,  however,  the  dominant  error  component 

is  due  to  the  occurrence  of  outliers  given  by  the  first  term  of  Equation  (67).  The  overall 

RMS  error  e  is  shown  in  Figure  9  as  a  fimction  of  SNR  with  the  sampling  rate  /,  =  10  KHz  ' 

and  the  value  of  m  equals  to  8.  Comparing  Figures  7  and  9,  the  effect  due  to  the  MSE  of 

the  DFT  approximation  can  be  easily  observed.  As  stated  earlier,  this  error  increasingly  t 

becomes  the  dominant  factor  as  the  SNR  increases.  This  is  because  the  variance  of  the 

ML  estimation  given  by  the  C-R  bound  is  inversely  proportional  to  the  SNR,  whereas  the 

DFT  iq>proximation  error  is  independent  of  SNR.  Hence,  with  increasing  values  of  SNR, 

the  MSE  of  the  ML  estimation  becomes  negligible  compared  to  the  MSE  of  the  DFT 

approximation  which  remains  basically  constant  for  a  given  value  of  m  and  N. 
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6.0  SUMMARY  AND  CONCLUSION 


In  this  paper,  we  have  studied  the  problem  of  estimating  the  signal  parameters, 
more  specifically  the  doppler  frequency,  from  a  finite  number  of  noisy  I  (in-phase)  and 
Q  (quadrature)  samples.  The  received  signal  is  2issumed  to  be  coherent  and  therefore, 
samples  ban  be  taken  at  regular  intervals.  In  addition,  the  signal  is  also  assumed  to  be 
corrupted  by  additive  white  Gaussian  noise.  The  equations  that  describe  the  Cramer-Rao 
bounds  which  give  the  minimum  variance  of  an  unbiased  estimator  are  derived.  In  the  case 
of  the  frequency  estimation,  the  C-R  lower  bound  is  shown  to  be  inversely  proportional  to 
SNR,  the  square  of  the  sampling  time-interval  T,  and  the  cube  of  the  number  of  samples 
N  observed.  As  for  the  phase  estimation,  it  is  found  that  the  C-R  bound  is  dependent  on 
the  time  fo  at  which  the  first  sample  is  taken.  The  lower  variance  bound  is  2dso  derived 
for  a  more  generalized  chirp  signal  whose  phase  is  modeled  2is  a  polynomial  of  order  p, 
where  p  >  1.  The  example  of  the  C-R  bounds  for  the  coefficients’  estimates  of  a  cubic 
polynomial  (p  =  3)  was  described  in  detail. 

The  maximum-likelihood  estimation  for  the  unknown  parameters  of  both  the 
simple  sinusoid  and  the  chirp  signal  are  derived.  The  relationship  of  the  ML  frequency 
estimate  to  the  discrete  Fourier  transform  is  then  shown.  The  ML  estimator  is  the  min¬ 
imum  variance  estimator  but,  however,  its  analysis  reveals  that  at  low  SNR,  it  is  also 
plagued  by  the  threshold  effect.  That  is,  there  is  a  range  of  low  SNR  for  which  the  MSE 
starts  increasing  very  rapidly  as  the  SNR  decreases.  This  is  because  of  the  large  error 
committed  by  the  ML  estimator  in  the  presence  of  an  outlier.  The  threshold  behaviour 
of  the  ML  estimation  was  discussed  in  Section  4.3. 

The  algorithm  for  the  implementation  of  the  ML  frequency  estimator  was  then 
considered.  It  consists  of  two  parts;  a  coarse  search  and  a  fine  search.  The  coarse  search 
3rields  a  rough  estimate  of  where  the  global  maximum  of  the  spectrum  is  located  and  this 
is  performed  using  the  fast  Fourier  transform  algorithm.  The  fine  search  uses  the  coarse 
estimate  as  the  initial  guess  in  order  to  give  a  more  accurate  result  of  the  frequency 
estimation.  It  can  be  performed  using  numerical  procedures  such  as  the  bisection  and  the 
secant  methods.  Errors  due  to  the  direct  estimation  of  the  frequency  parameter  by  the 
FFT  algorithm  was  also  analysed.  It  was  found  that  at  moderate  to  high  SNR  (above 
the  threshold  effect),  the  dominant  factor  in  the  overall  MSE  of  the  frequency  estimation 
is  the  DPT  approximation  to  the  ML  estimate.  The  accuracy  of  this  approximation  can 
be  increased  by  increasing  the  number  of  DFT  points  but  at  the  expense  of  more  number 
of  computations.  Once  the  number  of  DFT  points  is  chosen,  the  MSE  due  to  the  DFT 
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approximation  is  constant  and  independent  of  the  input  SNR. 

We  have  determined  the  minimum  variance  frequency  estimation  using  the  ML 
estimator.  From  the  expression  of  its  C-R  bound,  it  is  seen  that  this  estimation  can 
be  improved  by  increasing  the  effective  time  duration,  NT^,  of  the  received  signal.  The 
longer  the  time  duration,  the  more  accurate  the  frequency  of  the  signal  can  be  estimated. 
Using  the  proposed  method  for  estimating  the  frequency  of  the  signal  during  a  single 
pulse  of,  say,  one  microsecond,  with  a  sampling  rate  of  10*  Hz,  and  a  SNR  =  2Q dB,  the 
best  RMS  error  that  can  be  achieved  is  3.899  KHz.  Since  the  doppler  frequency  expected 
in  practice  is  of  the  order  of  a  few  hundred  hertz  or  less,  it  is  obvious  that  this  method 
cannot  measure  or  even  detect  any  doppler  shift  from  a  single  received  pulse.  The  large 
value  of  the  estimation  error  is  due  to  the  small  pulse  time  duration.  It  is  to  be  noted 
that  for  a  given  time  duration  r,  where  t  =  NT,,  the  accuracy  of  the  minimum  variance 
estimation  can  be  increased  by  increasing  N,  the  number  of  samples,  for  the  same  value 
of  T.  This  is  because,  from  the  expression  of  the  C-R  bound  in  (20),  the  MSE  is  inversely 
proportional  to  the  cube  of  N  (N^),  whereas  it  is  only  inversely  proportional  to  the  square 
of  the  sampling  interval  T,  (Tj*).  This  reduction  in  MSE  for  the  ML  estimation,  however, 
will  be  apparent  only  for  high  SNR  (above  the  threshold  effect)  zmd  it  will  be  of  almost 
no  effect  on  the  overall  MSE  if  the  DFT  approximation  is  used. 

6.1  Pulsed  Signals 

Finally,  some  examples  will  be  considered,  where  the  target  emitter  is  trans¬ 
mitting  either  a  burst  of  coherent  or  noncoherent  pulses  at  the  receiver.  The  goal  is  to 
determine  the  best  possible  estimation  that  can  be  achieved  with  a  finite  number  of  pulses. 
It  is  assumed  that  the  received  pulses  have  a  moderate  to  high  SNR  so  that  the  maximum 
likelihood  estimator  is  operating  above  the  threshold  value. 

In  the  case  of  the  coherent  pulses  where  basically  a  continuous  wave  signal  is 
being  modulated  by  a  square  waveform,  suppose  that  P  number  of  pulses,  each  with 
a  time  duration  of  r  seconds,  are  received.  The  pulsed  signal  is  sampled  at  a  rate  of 
f,  =  l/TfHz  which  implies  that  the  number  of  samples  per  pulse  is  N  =  tIT,  =  r/,. 
The  sampling  of  the  coherent  signal  can  be  performed  continuously  from  pulse  to  pulse 
such  that  the  integration  over  the  P  pulses  is  equivalent  to  the  sampling  of  a  CW  signal 
of  time-length  Pt  seconds.  The  MSE  of  the  ML  angulau:  frequency  estimation  is  then 
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obtained  from 


6  ^  6 
{SNR)TjNPiN^P^  -  1)  {SNR)T^N^P^ 


for  TV,  P  1. 


(69) 


Comparing  the  above  equation  with  that  of  (20),  it  is  seen  that  by  having  P  number  of 
coherent  pulses,  the  MSE  is  improved  by  a  factor  of  P^.  Thus,  in  terms  of  the  RMS 
error,  the  improvement  factor  is  of  the  order  of  P^^^.  For  the  example  used  earlier,  where 
r  =  1/is,  /,  =s  l(fi  Hz,  and  SNR  =  20  dB,  the  incretise  in  accuracy  of  the  frequency 
estimation  when  P  =  10  is  by  a  factor  of  31.62  that  yields  an  RMS  error  of  123.28 /fe. 
Similarly,  if  we  require  the  frequency  estimate  to  have  an  RMS  error  of  less  than  10  Hz, 
then  for  the  same  values  of  r,  /„  and  SNR,  the  number  of  coherent  pulses  P  needed  to 
achieve  this  accuracy  is  P  >  54. 


For  noncoherent  pulses,  where  every  received  pulse  is  independent  of  each  other, 
we  cannot  combine  the  pulses  together  to  form  a  CW  wave  as  in  the  coherent  case.  We 
have  to  treat  each  pulse  separately  and,  therefore,  the  best  we  can  do  in  order  to  improve 
the  accuracy  of  the  frequency  estimation  is  to  have  an  average  estimate  of  all  the  received 
pulses.  That  is, 

^  1=1 

where  /,•  is  the  frequency  estimate  of  the  ith  pulse.  The  variance  of  the  average  estimate 
/,  is  then  reduced  by  a  factor  of  P,  the  number  of  pulses,  i.e., 

Var(/)  =  -  _  i)  C^O) 

which  implies  that  the  RMS  error  of  the  frequency  estimate  is  decreased  by  P^^*.  Now,  in 
order  to  have  an  accuracy  of  10  Hz  or  less,  with  each  pulse  having  the  values  of  r  =  1  /is, 
SNR  =  20  dB,  and  /,  =  l(PHz,  we  need  P  >  151,997  pulses  which  is  too  large  for 
practical  purposes.  Note  that  an  accuracy  of  the  order  of  10  Hz  is  chosen  because  that  is 
basically  what  is  required  in  practice  for  doppler  detection  and  measurement. 

One  way  of  reducing  the  number  of  pulses  P  for  the  same  accuracy  is  to  have 
more  samples  per  pulse,  i.e.,  increase  the  value  of  N  by  sampling  at  a  faster  rate,  as  was 
explained  before.  If  the  number  of  samples  is  augmented  from  N  to  QN,  where  Q  >  1, 
then  the  MSE  is  approximately  reduced  by  Q  and  the  RMS  error  by  As  for  the 

number  of  noncoherent  pulses  needed,  P  is  decreased  to  P/Q  and  that  of  coherent  pulses  to 
P/Q^^^  for  the  same  MSE  or  RMS  error  value.  In  the  example  of  the  noncoherent  pulses, 
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if  the  sampling  rate  is  increased  from  10*  to  10® //s  for  the  same  value  of  t,  {Q  =  10), 
then  for  the  same  accuracy  of  iOHz  or  less,  at  least  15,199  pulses  are  required,  which  is 
still  very  large.  Comparing  the  noncoherent  pulses  to  the  case  of  the  coherent  ones,  where 
the  number  of  pulses  required  is  in  the  order  of  hundred  thousand  instead  of  only  a  few 
tens,  it  can  be  deduced  that  it  is  only  practically  feasible  for  an  ESM  system  to  measure 
the  change  in  doppler  frequency  from  the  coherent  pulses  and  not  from  the  noncoherent 
pulses,  since  in  general,  the  number  of  pulses  received  during  a  burst  is  usually  a  few 
hundreds  or  less. 

In  the  case  where  the  demodulated  pulse  signal’s  frequency,  /  is  less  than  the 
reciprocal  of  the  pulse  duration  r  and  is  less  than  the  pulse  repetition  frequency  (PRF), 
i.e.,  /  <  1/t  and  /  <  PRF,  each  pulse  can  be  taken  as  a  discrete  sample  with  the 
PRF  as  the  sampling  rate.  However,  this  can  be  done  only  if  the  pulse  signal  is  coherent. 
For  a  PRF  of  10  KHz  and  a  SNR  of  20/,  dP,  using  the  expression  of  the  C-R  bound 
in  (20),  an  accuracy  of  10  Hz  or  less  can  be  achieved  with  at  least  12  pulses.  For  the 
frequency  /  <  1/t  but  greater  than  the  PRF  (/  >  PRF),  each  pulse  cannot  be  sampled 
only  once,  because  this  will  lead  to  ambiguities  in  determining  the  true  frequency.  More 
specifically,  there  will  be  a  difference  of  n(PRF)  between  the  true  frequency  and  the 
estimated  frequency,  where  n  =  1,  2,  3  •  •  • .  The  three  different  frequency  ranges  of  the 
pulse  signal  from  the  output  of  the  demodulator  are  shown  in  Figure  10.  For  /  >  1/t, 
the  signal  is  readily  discerned  from  the  information  contained  in  a  single  pulse,  whereas, 
for  /  <  1/r,  the  pulses  are  modulated  with  the  amplitude  of  the  signal  waveform. 
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Figure  10;  (a)  pulse  train  for  signal  frequency  /  >  7;  (b)  pulse  train  for  signal  frequency 
/  <  7  and  /  <  PRF]  (c)  pulse  train  for  signal  frequency  /d  <  7  and 
/d  >  PBF. 
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